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Abstract: Modified boundary conditions (MBC) and a multilayer 

approach (MA) are proposed as fast and efficient numerical methods for 
the design of ID photonic structures with rough interfaces. These methods 
are applicable for the structures, composed of materials with arbitrary 
permittivity tensor. MBC and MA are numerically validated on different 
types of interface roughness and permittivities of the constituent materials. 
The proposed methods can be combined with the 4x4 scattering matrix 
method as a field solver and an evolutionary strategy as an optimizer. The 
resulted optimization procedure is fast, accurate, numerically stable and can 
be used to design structures for various applications. 
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1. Introduction 

Photonic crystals are of high interest for many practical applications [1, 2, 3]. One-dimensional 
photonic crystals are known for a long time and have been studied thoroughly [4, 5, 1, 6]. A 
lot of effort was also put into the optimization of finite multilayer films [7, 8, 9, 10]. Different 
optimization procedures, such as the quasi-Newton algorithm [8] or genetic algorithms [9, 10] 
have been applied to improve the performance of finite multilayer mirrors. 

Some applications require the optimization of the structure for high reflection or anti- 
reflection in more than one frequency band. Reentry [11] is an example of such an application, 
where reflection in two ranges is desired. 

In practice roughness of the interfaces is inevitable and causes severe numerical problems, 
because it destroys the symmetry of the structure. In order to avoid time-consuming 3D simula- 
tions, we propose a multilayer approach (MA) or modified boundary conditions (MBC) to take 
roughness into account. These approaches are combined with the 4x4 scattering matrix method 
as field solver. The resulting procedure is fast, numerically stable, and easy to implement. De- 
tails of the scattering matrix formalism can be found in [12, 13, 14], 

This paper is organized as follows. In Section 2 we illustrate the optimization procedure on 
a dual band mirror with flat interfaces. In Section 3 and Section 4 we describe MA and MBC 
for rough interfaces in detail. In Section 5 we provide numerical results of MBC and MA for 
rough interfaces. A sinusoidal and a random hilly profiles of roughness are taken as examples. 
Results for isotropic and anisotropic uniaxial crystals are shown. 

2. Optimal design of a dual band omnidirectional mirror 

In most cases it is not possible to find an optimal design of a structure analytically. Various 
optimizers [15, 16] can be helpful to find good solutions numerically. 

If a layered structure consists of N layers, one has (at least) N real-valued parameters (thick- 
ness of each layer) to be optimized. The search space should be defined for each parameter 
within a reasonable interval [c/ m j n , aWxL For such optimizations, various algorithms, such as 
genetic algorithms (GA), micro-genetic algorithms (MGA), or evolutionary strategies (ES) may 
be applied. The optimization task becomes difficult not only when N is high, but also when the 
fitness function has a complicated behavior, e.g. when it has many local optima. According to 
our experience [17] ES is very powerful in real parameter optimization problems and outper- 
forms genetic algorithm (GA), particle swarm optimization (PSO) etc. in most cases. Therefore, 
we applied it to the multilayer problem without extensive testing. We used an ( m+n ) ES with 
adaptive mutation for the optimization in the following example. Here m is the initial number 


of parents, n is the number of children created in each generation. In the following example we 
used m = 5 and n = 40. 

To illustrate the optimization procedure, we optimize the omnidirectional reflection of the 
structure in two bands. Hereafter, we denote the reflection coefficients: R ss , R sp , R ps , and R pp . 
The first and second indices correspond to polarization of the incident and reflected light re- 
spectively. The s- and p-polarized light correspond to electric or magnetic field parallel to layers 
respectively. The multilayered structure to be optimized consists of N= 14 layers, embedded be- 
tween air interfaces. It is composed of two alternating materials with permittivity tensors £i and 
£ 2 : 
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We take these permittivities just as examples. Any other materials with full permittivity ten- 
sors could be used instead. In this example, we consider flat interfaces between the layers. In 
case of rough interfaces MBC or MA can be used (see Section 3 and Section 4). MBC and MA 
are applicable when the characteristic size of the roughness is much smaller than the wavelength 
in the material (5 << A). 

As a necessary step for the optimization, we must define the fitness function, i.e. the goal 
of the optimization. We optimize the omnidirectional reflection of the unpolarized light in two 
bands B\ := 800 — 1300 meV and Ih '■= 2000 — 2500 meV. Therefore, we define the fitness 
function as the averaged reflection: 

Fitness = (R(di,d 2l ... ,dff)) = ^ ^ ^ ^ 1 • (1) 

f f 0.5 (R ss FR S p ~hRps FR pp )d cod 6 

(R U 2(d h d 2 ,... ,d N )) = — J jd(0de <1- (2) 
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where the reflection coefficients are integrated over the angles of incidence 0 £ [0, and 
energy bands B\ and Ih- The thickness of the i-th layer (i = 1 . .... A) is t/,. We now need to 
maximize the fitness function Eq. (1) of N variables. 


(a) Reflectivity=0.5(Rss+Rsp+Rpp+Rps) (b) 



Fig. 1. (a) Dependence of the total reflectivity of the unpolarized light on the energy of 
photons and the angle of incidence 8 . (b) Thicknesses of individual layers of the optimal 
structure. Total number of layers N=14. 


The results of ES the optimization are shown in Fig. 1 . One may see that a high reflectivity is 
achieved in the B\ and B 2 bands. Figure l.b shows the geometric profile of the obtained optimal 
structure. 



In the optimization procedure we limited the number of fitness calls to 14200. As a conse- 
quence it is not guaranteed that the global optimum is found. As long as a sufficiently good 
solution is found, this is acceptable in practice. The required optimization time was several 
hours in Matlab on Intel(R) Core(TM) i7 CPU 2.67GHz with 8 GB of RAM. 

The interfaces of fabricated multilayered films are not ideally flat. Roughness between inter- 
faces can be a result of fabrication inaccuracies or interdiffusion, if the structure operates at high 
temperatures. Therefore, it is desired to take roughness into account at least approximately. We 
propose two ways for doing this: 1 ) multilayer approach (MA) 2) modified boundary conditions 
(MBC). These methods work very fast, this is essential for the optimization problems. 


3. Multilayer approach for rough boundaries 

Model 

Geometry. The rough boundary is located between two subspaces: z < £ = £i, ju = 1, and 

Z > +5, £ = £ 2 , jti = 1. The region |z| < 8 is occupied with the rough boundary, see Fig. 2. Both 
sizes of the roughness 28 and / are supposed to be small in comparison with the wavelength. 
The volume fractions of two components are j\ = V\ jV and jo = V 2 /V, where V is the total 
volume of the rough region. 



Fig. 2. Rough boundary between interfaces with permittivities £\ and £ 2 . 

We can divide the region — 8 < z < 8 into M layers and assign an effective permittivity G, ir 
(i=l,...,M ) to each of these layers. However, standard formulas for effective permittivity, such 
as Maxwell-Garnett orBruggeman, are not valid in the general case of arbitrary anisotropic ma- 
terials and geometric profile of roughness. Therefore, we can only use the bounds for effective 
permittivity: 


(' V '■ — — Emin < £eff < £max — (/l £l + fl&l) (3) 

V £ 1 e 2 / 

These bounds for effective permittivity were first derived by Wiener [18] and they are valid for 
arbitrary anisotropic materials and shapes of the inclusions. 

As an effective permittivity for each layer, the mean value of the two bounds can be taken: 
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4. Modified boundary conditions for rough boundaries 

Model 

Geometry. The rough boundary is located between two subspaces: z < — 8 , £ = £i, jJL = 1, and 
Z > +5, £ = £ 2 , H = 1. The region |z| < 8 is occupied with the rough boundary, see Fig. 2. We 


suppose that the dielectric permittivity is a function of z in this region. Thus, 


£ 1 , 

if z < - 8 

e(z) 

if z| < 8 

£ 2 , 

if z>+5. 


The wave is incident from z = — and propagates in the XZ plane, thus (d /dy = 0). The angle 
of light incidence is 6 with respect to Z axis. 

Let us choose the rectangular / x 28 contour at the rough boundary. Then the Maxwell’s 
equations for the circulation around this contours are: 
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We can write the left part of the first equation as: 
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because B„ is continuous at the boundary and can be considered as constant there. Thus we 
obtain: 

fndl = [B^\-8)--a^(5)]l+[B n (l/2)-R n (-l/2)]J S s 

The right part of the first equation (5) can be written as: 
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Assuming 8 <C A and l <sjC A, we obtain 
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Here all fields are written at z = 0. The first and second equations originate from the first 
equation (5) being written for the different components of the magnetic field. The third and 
forth equations are derived from the second equation of (5). 


The boundary conditions Eq. (6) can be simplified, if we consider a wave of a certain polar- 
ization. For E-waves we have 


E = (0, E y , 0) 

H = (H x , 0, H z ). 

The second and third equations are satisfied identically, while the first and the last equations 
reduce to 


£y (1) -£y (2) = 0 (7) 

dE™ dE y 2) _ r (i) 
dz dz c 2 1 y 


Similarly, we have for H-waves 


E = (E x , 0, E z ) 
H= (0, H y , 0) 


and we obtain 


Here 
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The generalization of the boundary conditions Eq. (7) and Eq. (8) for the anisotropic struc- 
tures leads to very lengthy expressions. However, they are essentially simplified for uniaxial 
crystals, where the dielectric permittivity tensor is diagonal and £n = £22 ^ £33. Then bound- 
ary conditions Eq. (7) and Eq. (8) hold with £1 = e{^, £2 = eJ/* and 


h =5(^£ i ( ] ) +£ 1 (2) ) - j e n {z)dz 

f s dz 8 8 
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Here the superscripts (1) and (2) correspond to the first and second bordered material, respec- 
tively. 

The second order correction on 5/A needs an additional expansion of the fields in the inte- 
grands of Eq. (5). This can be done also by integration of the wave equations over the interface 
region |z| < 8, see Ref. [19]. The wave equation for the E-wave in this region is: 


d 2 Ey d 2 Ey 

dz 2 + dx 2 


ft) 


+ £(z)— Y E y ~ 0- 


(9) 


The integration of this equation over the 5-region at the boundary (see Ref. [19] for the details 
of calculation) leads to the following boundary conditions for the electric fields E\ 2 and their 
derivatives E[ 2 = dE\ 2 /dz at each side of the boundary (we omit the subscript y for reasons 
of simplicity): 



Here a and are the parameters, which are independent of the wavelength and incidence angle. 
The wave equation for the H-wave in 5 region is 
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if £33 is independent of x, then Eq. (11) can be written as: 
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Expansion similar to Eq. (6), but to the second order on 5/A, leads to a = t\, = (£2 — £1 )5 2 /2 

and y = f 2 - Here the parameters a and /j are the same as in Eq. (10). 


5. Numerical results of modified boundary conditions and multilayer approach 

We demonstrate the application of MBC and MA for rough interfaces by the following exam- 
ples. 

5.1. Sinusoidal roughness profile 

Consider the reflection from a slab with a sinusoidal grating on the top, see Fig. 3. a. The sinu- 
soidal grating emulates a rough boundary between air and slab. The permittivities of the slab 
and of the sinusoidal grating are the same £ = 5. The structure is embedded in air. The dielectric 
permittivity in the region —5/2 < z < 5/2 is: 

f *2=5, ifz>5(±— sin(f)sin(f)^ 

\ *1 = 1, ifz<5(£-sin(f )sin(S )j 

Period and height of the grating are much smaller than the considered wavelength. As a conse- 
quence, the dependence of the reflectivity on the azimuthal angle cp is almost negligible, despite 
the absence of the symmetry with respect to cp. 

We calculate the reflection spectra of the structure in four different ways. The results are 
shown in Fig. 4. a and Fig. 4.b. The red curve represents the solution obtained with CST Studio 
(Frequency Domain Solver). We refer to it as the ’’reference solution”. The black dash-dotted 
curve represents the reflection from a flat slab with thickness L + h/2 = 105 nm. One can 
see that the reflection from a flat interface deviates significantly from the ’’reference solution”. 



Fig. 3. (a) Sinusoidal grating on top of the slab with thickness L=100 nm. Height and 
periods of the grating are h = 10 nm, d x = d y = 10 nm. (b) Random hilly roughness on top 
of the slab with thickness L=100 nm. Periods of the unit cell are d x = d y = 50 nm, maximum 
height of the hills is /; max = 10 nm. (c) Multilayered structure with rough interfaces 1-3. 


This is especially pronounced for higher energies. The blue curve corresponds to the multilayer 
approach, with M = 6 layers in the sinusoidal region —h/2 <z< h/2. The green curve corre- 
sponds to the MBC solution (Eq. (10) and Eq. (12)), with optimal parameters a [cm], /j [cm 2 ] 
and y [cm] . 

In Fig. 4 one can see, that MA is in very good agreement with the reference solution for both 
polarizations and for all angles of incidence 0. MBC also is in good agreement. The parameters 
a, ji and 7 of the MBC are independent of the frequency and incidence angle. These parameters 
can be found either by fitting to the reference solution, as done in this paper, or by fitting to 
experimental reflection/transmission spectra of multilayers. 

Now we consider the reflection from the same structure, but composed of anisotropic material 
with £n = £22 = 7, £33 = 9. It should be noted that the application of MBC is very much 
simplified in the case of isotropic, or uniaxial crystals (£n = £22 ^ £ 33 ). While MA is equally 
easy to apply for arbitrary permittivity tensors. 

The results of MA and MBC are shown in Fig. 4.c and Fig. 4.d. It can be seen that MA (blue 
curve) is in very good agreement with the reference solution (red curve). MBC also show a very 
good agreement. Thus the application of MA and MBC to the considered sinusoidal roughness 
profile is justified. 

5.2. Random hilly roughness 

In order to demonstrate that MA and MBC methods are not restricted to sinusoidal approxima- 
tion of rough interface, we consider the reflection from the slab with random hills (polygons) 
on top, see Fig. 3.b. Maximum height of the hills does not exceed /i max = 10 nm. The roughness 
profile is periodic, with periods d x = d y = 50 nm. The characteristic size of the roughness is 
much smaller than the wavelength, therefore the dependence of the reflectivity on the azimuthal 
angle <p is practically negligible. 

As in Sec. 5.1 for the sinusoidal roughness, we consider separately isotropic and anisotropic 
materials. The results of MA and MBC for an isotropic material with £ = 5 are shown in 
Fig. 5. a and Fig. 5.b. The results for an uniaxial crystal with S\ 1 = £22 = 7, £33 = 9 are shown in 
Fig. 5.c and Fig. 5.d. Evidently, a good agreement of MA and MBC with the reference solution 
is observed. 

If a structure consists of N layers with rough interfaces, then either MA or MBC can be 



Fig. 4. Reflectivity of the structure shown in Fig. 3. a. (a) and (b) correspond to an isotropic 
material with e = 5, (c) and (d) correspond to an anisotropic material with £n = £22 = 
7, £33 = 9. Red curve - reference solution (CST Studio), black dash-dotted curve - flat slab 
with thickness 105 nm. blue curve - multilayer approach (MA). green curve - MBC between 
air and slab. 


applied to each of the N-l interfaces. This is demonstrated in the next example. 

5.3. Multilayers with rough interfaces 

In this final example, we apply MA and MBC to a multilayered structure with rough interfaces. 
The structure is shown in Fig. 3.c. Structure is composed of two materials with permittivities 
£i= 2 and £2 = 7, high index material is on top. Interfaces 1-3 are rough, with random hilly 
roughness profiles, as in Fig. 3.b. Roughness of each interface was generated randomly. The 
maximum height of the hills does not exceed /z max = 10 nm. The structure is periodic with 
periods d x = d y = 50 nm. 

Reflection spectra of a normally incident light is shown in Fig. 6. The red curve corresponds 
to the ’’reference solution”, obtained with the CST Studio. Some resonances are observed for 
energies above 9 eV. The reason is that the size of roughness becomes comparable with the 
wavelength. The black dash-dotted curve corresponds to flat interfaces. The blue curve corre- 
sponds to MA with M = 6 layers in the interface regions. The green curve on the plot corre- 
sponds to MBC at each rough interface. One can see that MA and MBC in general demonstrate 
a better agreement with ’’reference solution”. This is especially pronounced in two regions, 
which are marked with arrows in the Fig. 6. 

It should be noted that MA and MBC combined with scattering matrix method work very 
fast. Calculation of one frequency point of the structure shown in Fig. 3.c requires only about 
1 ms with the scattering matrix code and about 1 min with the CST Studio (Frequency Domain 
Solver). Calculations were performed on Intel(R) Core(TM) i7 CPU 2.67GHz with 8 GB of 
RAM. 



Fig. 5. Reflectivity of the structure shown in Fig. 3.b. (a) and (b) correspond to an isotropic 
material with e = 5, (c) and (d) correspond to an anisotropic material with £n = £22 = 
7, £33 = 9. Red curve - reference solution (CST Studio), black dash-dotted curve - flat slab 
with thickness 105 nm. blue curve - multilayer approach (MA). green curve - MBC between 
air and slab. 



Fig. 6 . Reflectivity of the structure shown in Fig. 3.c for normally incident light. Permit- 
tivities of materials are £1 = 2 and £ 2=7 (on top). Red curve - reference solution (CST 
Studio), black dash-dotted curve - flat interfaces, blue curve - multilayer approach (MA) at 
interfaces, green curve - MBC at interfaces. 


6. Conclusions 

We presented modified boundary conditions (MBC) and a multilayer approach (MA) as fast and 
efficient numerical methods for the design of layered structures with rough interfaces. These 
methods are applicable for arbitrary anisotropic materials. MA and MBC provide a good agree- 
ment with the reference solutions for sinusoidal and random hilly roughness profiles. Combined 
with the scattering matrix technique as field solver and an evolutionary strategy as numerical 
optimizer these methods reveal a fast way of optimizing ID photonic structures for various 


applications. 


